/*
   FieldVis - A Tool for Visualizing Astrophysical Magnetohydrodynamic Flow
   
   Copyright (C) 2006 - 2008     FieldVis Development Team
  
   This program is free software; you can redistribute it and/or modify
   it under the terms of the GNU General Public License as published by
   the Free Software Foundation; either version 2 of the License, or
   (at your option) any later version.
  
   This program is distributed in the hope that it will be useful,
   but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   GNU General Public License for more details.
  
   You should have received a copy of the GNU General Public License along
   with this program; if not, write to the Free Software Foundation, Inc.,
   51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
 
   File     : StreamlinePath.cpp
   Author   : Blayne Field
   
   Purpose  : Implementation file for the StreamlinePath class.  This class
              holds the coordinates and directions of a streamline.
*/

#include "StreamlinePath.h"


StreamlinePath::StreamlinePath()
   : coords(NULL), directions(NULL), pointsCalculated(0),
     startPointIndex(0), exitsGrid(false)
{
}

StreamlinePath::StreamlinePath(const StreamlinePath &slp)
{
   pointsCalculated = slp.pointsCalculated;
   startPointIndex = slp.startPointIndex;
   exitsGrid = slp.exitsGrid;

   coords = new Point[pointsCalculated];
   directions = new Vector3D[pointsCalculated];

   for (int i = 0; i < pointsCalculated; i++)
   {
      coords[i] = slp.coords[i];
      directions[i] = slp.directions[i];
   }
}

// stitches together two paths originating from the same point
// in opposite directions
StreamlinePath::StreamlinePath(StreamlinePath forward, StreamlinePath backward)
{
   // figure out how many points there are
   int forwardPoints = forward.pointsCalculated;
   int backwardPoints = backward.pointsCalculated;
   pointsCalculated = forwardPoints + backwardPoints - 1;
   coords = new Point[pointsCalculated];
   directions = new Vector3D[pointsCalculated];

   int newPathIndex = 0;

   // start by adding the points from the backward streamline
   for (int i = backwardPoints - 1; i >= 1; i--)
   {
      coords[newPathIndex] = backward.coords[i];
      directions[newPathIndex] = backward.directions[i-1] * -1;
      newPathIndex++;
   }

   // add the points from the forward streamline
   // (including the common start point)
   for (int i = 0; i < forwardPoints; i++)
   {
      coords[newPathIndex] = forward.coords[i];
      directions[newPathIndex] = forward.directions[i];
      newPathIndex++;
   }  

   startPointIndex = backwardPoints;
   exitsGrid = forward.exitsGrid;
}

StreamlinePath::StreamlinePath(const std::list<BoundVector3D> &locs, const bool exits)
   : startPointIndex(0), exitsGrid(exits)
{
   pointsCalculated = locs.size();

   coords = new Point[pointsCalculated];
   directions = new Vector3D[pointsCalculated];

   std::list<BoundVector3D>::const_iterator iter = locs.begin();

   for (int i = 0; iter != locs.end(); iter++, i++)
   {
      coords[i] = iter->base;
      directions[i] = iter->direction;
   }
}

StreamlinePath::~StreamlinePath()
{
   if (pointsCalculated > 0)
   {
      delete [] coords;
      delete [] directions;
      pointsCalculated = 0;
   }
}


void StreamlinePath::ScaleByFactor(float factor)
{
   for (int i = 0; i < pointsCalculated; i++)
   {
      coords[i] = coords[i] * factor;
   }
}

void StreamlinePath::SmoothPath(VectorField *field)
{
   int oldStartIndex = startPointIndex;

   if (pointsCalculated == 1)
      return;

   CubicSplinePolynomial3D *splines =
   new CubicSplinePolynomial3D[pointsCalculated - 1];
   int *pointsToCompute = new int[pointsCalculated - 1];
   int pointCount = 0;
   float stepSize = 1.0 / field->GetScaleFactor();
   float curvature;  

   // figure out how many new points are going to be needed
   // and calculate the splines as you go
   float scaleFactor = field->GetScaleFactor();
   for (int i = 0; i < pointsCalculated - 1; i++)
   {
      // scale the points
      Point p0 = coords[i] * scaleFactor;
      Point p1 = coords[i+1] * scaleFactor;

      // get the derivatives
      Vector3D d0 = field->GetInterpolatedVector(p0);
      Vector3D d1 = field->GetInterpolatedVector(p1);

      // make the spline
      splines[i] = CubicSplinePolynomial3D(p0, p1, d0, d1);
      curvature = splines[i].AverageCurvature();

      int distSteps = (int) (Point::DistanceBetween(coords[i], coords[i+1]) / stepSize);
      int maxSteps = 20 * distSteps + 1;
      int curvSteps = (int) (curvature * 10) + 1;
      pointsToCompute[i] =  (distSteps + 1 > curvSteps ? distSteps + 1:
                             (curvSteps < maxSteps ? curvSteps : maxSteps));    

      pointCount += pointsToCompute[i];
      if (i < oldStartIndex)
         startPointIndex = pointCount;
   }

   // get a new chunk of memory
   if (coords)
      delete [] coords;
   if (directions)
      delete [] directions;
   coords = new Point[pointCount];
   directions = new Vector3D[pointCount];

   int newPointIndex = 0;
   // interpolate each step
   for (int i = 0; i < pointsCalculated - 1; i++)
   {
      for (int j = 0; j < pointsToCompute[i]; j++)
      {
         coords[newPointIndex] = splines[i].EvaluateAt((float)j / pointsToCompute[i]);

         if (newPointIndex > 0)
            directions[newPointIndex - 1] =
            Point::VectorFromTo(coords[newPointIndex-1], coords[newPointIndex]);

         newPointIndex++;
      }
   }
   pointsCalculated = newPointIndex;
   ScaleByFactor(1 / scaleFactor);
   directions[newPointIndex - 1] = directions[newPointIndex - 2];

   delete [] splines;
   delete [] pointsToCompute;
}

StreamlinePath & StreamlinePath::operator=(const StreamlinePath & toCopy)
{ 
   // copy the object if it is different from this one
   if (this != &toCopy)
   {
      if (pointsCalculated > 0)
      {
         if (coords)
            delete [] coords;
         if (directions)
            delete [] directions;
      }

      pointsCalculated = toCopy.pointsCalculated;
      startPointIndex = toCopy.startPointIndex;
      exitsGrid = toCopy.exitsGrid;

      coords = new Point[pointsCalculated];
      directions = new Vector3D[pointsCalculated];

      for (int i = 0; i < pointsCalculated; i++)
      {
         coords[i] = toCopy.coords[i];
         directions[i] = toCopy.directions[i];
      }
   }

   return *this;
}

